Climate-smart forestry through innovative wood products and commercial afforestation and reforestation on marginal land

Significance Afforestation and reforestation (AR) are nature-based solutions to climate change. However, the greenhouse gas (GHG) mitigation efficacy of protection or commercial AR is under debate. This study develops a dynamic life cycle assessment to quantify the GHG mitigation potential of protection and commercial AR on marginal land in the southeastern United States. We found that commercial AR with cross-laminated timber and biochar production generally mitigates more GHGs across 100 y than protection AR and commercial AR with traditional lumber production. Protection AR could mitigates more GHGs in a shorter timeframe (≤50 y). These results highlight the role of synergizing protection AR, innovative wood utilization, and strategic forest plantation management in supporting short- and long-term climate change mitigation goals.


SI Note 1: Summary of Methodology
This study presents a cradle-to-grave multi-scale dynamic LCA that integrates a GIS model, a forest growth model, SOC model, and process-based models of manufacturing wood products (i.e., lumber, CLT, and biochar) to quantify the carbon implications of afforestation and reforestation (AR) on the marginal land in the southeastern US over 100 years. The GIS model extracted and processed marginal land layers from a study (1) that identified AR opportunities across the US. The carbon sequestration by forest and biomass carbon was simulated by the forest growth and yield model over 100 years. Then carbon emissions from soil were simulated by the SOC model RothC (2) with the carbon input from the forest (e.g., litterfall, roots, residues, and snags). The GHG emissions by the production of three wood products were simulated by process-based models. This study includes the carbon associated with the end of life (EOL) of wood products as well as the upstream production of fuels, chemicals, and electricity consumed within the system boundary. Scenarios are established to explore the impacts of various forest plantation management strategies and product choices. The functional unit is the whole available marginal land for AR in the southeastern US (2.1 Mha). The Intergovernmental Panel on Climate Change (IPCC) AR6 GWP-100 factors are used to convert GHG emissions to CO2 equivalent basis (3). Detailed methods are discussed in the sections below (SI Notes 2 to 9).
In this cradle-to-grave LCA, the net GHG balances ( ) and the total carbon stock ( ) are calculated by equations (S1) and (S2), respectively. Fig. S8 shows the carbon flows associated with different life-cycle stages and carbon pools.
Where − 2 (estimated by forest growth model) is the accumulative CO2 sequestration by forest growth, (estimated by RothC model) stands for accumulative CO2 emissions from the forest that includes the biogenic emissions from biomass decomposition in soil and forest floor, represent fossil GHG emissions caused by fuels used in forest operation. and are biogenic and fossil GHG emissions at the stage of production (producing wood products and biochar), respectively. and are biogenic and fossil GHG emissions at the end-of-life stage of the products. In equation (S2), (estimated by forest growth model) is the carbon stock in forest standing, (estimated by RothC model) stands for carbon stock in soil and forest floor, and and represent carbon stock in products (wood products and biochar) and product end-of-life, respectively. All GHG emissions and carbon stock changes during wood product production and wood product end-of-life along with fossil GHG emissions by forest operation are estimated by the LCA model for wood products manufacturing.
For baseline (protection AR), the net GHG balance ( ) can be adapted shown as equation (S3) based on equation (S1) because there is no log harvesting for wood product manufacturing (or say , , , , , are zero).
Biogenic carbon mass balances of the forest system (I in Figure S8): inputs = outputs + carbon stock changes. Inputs are carbon uptake ( ), outputs are logs and residues and carbon emissions from forests ( ). The carbon stock changes include changes in carbon stock in standing forests ( ∆ ) and soil and forest floor (∆ ).
Based on equation (S4) and the assumption that all carbon losses from forest systems are emitted as CO2, we can derive the equation (S5) for baseline production AR.
& are intermediate flows that go to wood product manufacturing, which are zero for the baseline production AR. − 44 12 • �∆ According to equation (S5), the net differences between biogenic GHG emissions and CO2 sequestration from the forest system (− 2 + ) can be calculated using the carbon stock changes of standing forests and soil & floors (− 44 12 • �∆ + ∆ �. We use this method (S5) and (S3) to calculate the net GHG balances of the baseline protection AR because cannot be directly estimated for a 100 year time frame -the forest growth model has a maximum modelling period of 35 years, without the annual carbon input from year 36 to 100, RothC cannot simulate GHG emissions from the soil and forest floor. The calculation of carbon stock change in the forest standing ∆ are documented in S1 Note 2 and carbon stock change in soil and forest floor ∆ are documented in SI Note 3.

SI Note 2. Forest growth and yield
In this study, we used the 1996 Plantation Management Research Cooperative (PMRC) whole stand growth and yield model to simulate forest growth for site-prepared loblolly pine (Pinus taeda L.) plantations across the three prominent physiographic regions in the southern US: piedmont, upper coastal plain, and lower coastal plain in the southeastern US (4). The data for establishing the models come from hundreds of permeant plots in Alabama, Georgia, North Carolina, South Carolina, and Virginia, the first of which was established in 1977. Those plots were carefully established across a wide geographic range and re-measured throughout the rotation of the stand. Each individual tree in the permanent plot was measured and those were compiled into stand-level estimates of dominant height, basal area, and tree density. Those stand-level estimates were used to build the 1996 PMRC whole stand models used in our analysis. The model consists of equations to predict the whole stand: 1) dominant height/site index, 2) survival rate as trees per acre, 3) basal area per acre, 4) yield per acre, and 5) yield breakdown. Detailed equations can be found in the PMRC Technical Report (4). In addition to the whole stand estimation model, thinned plantation and the effects of mid-rotation fertilization with N and P are considered for the growth and yield simulation in the PMRC model. Thinned basal area is estimated as a function of the number of trees thinned. A thinned growth response in terms of per-acre basal area is formulated by comparing the basal area growth of a thinned plantation to that of an un-thinned counterpart of the same age, dominant height, and number of trees per acre. The effect of fertilization is estimated by accounting for a growth response in the dominant height and per-acre basal area growth equations. The growth response is estimated as a function of pounds of elemental nitrogen per acre, whether Phosphorus was applied, and the number of years since treatment.
The main input variable for the model under each scenario is site index base age 25 for loblolly pine with a range from 50 to 105 feet. The GIS data of site index for loblolly pine was extracted from the Gridded Soil Survey Geographic (gSSURGO) Database (5). Key variables that differ by scenario include: planning horizon (years), stand density (trees/acre), whether and when thinning is included, and density removed by thinning. Other common variables are whether fertilization is included, fertilization age, and fertilization application rate. The settings, input parameters, and variables for the model in this study are shown in Table S13 The annual total dry mass of the whole living stand, snags, and removed trees from thinning of loblolly pine for a specific site index was estimated using equations (S6-S8) from literature where the dry mass metric ton of foliage per ha ( , Mg/ha), branch ( , Mg/ha), stem ( , Mg/ha), and root (Mg/ha) can be calculated with modeled diameter outside-bark at breast height 1.37m (4.5 ft) ( ℎ, cm), total height ( , m), age ( ), diameter distribution, and annual standing density that is generated from PMRC 1996 model (7). Other outputs from PMRC model include basal area ( , m 2 /ha) and quadratic mean diameter ( , cm). A ℎ to ratio of 1.04 based on the literature (7) was used to calculate the ℎ.
Where represents living tree density or numbers of snags and thinned trees each year (number of trees per ha), 1000 is unit conversion factor from metric ton to kg, , , and (n is from 1 to 5 and m is from 1 to 4) are parameters used in equations (S6-S8) for loblolly pine and are listed in Table S14 (7).
The annual dry mass of litterfall biomass was calculated based on the needlefall ( , Mg/ha) to litterfall ( , Mg/ha) ratio ( / ) that was calculated by equation (S9) from literature (8). Annual day mass needlefall biomass was estimated by equations (S10-S15) from literature (8). Root mass (Mg/ha) is assumed to be 50% of stem mass, according to literature (9).
Where is leaf area index, is stand density index in metric unit, represents living stand density (number of trees per ha), is site index (m), , , 1 , 4 , and 1 (n is from 1 to 5) are parameters used in equations (S9-S15) for loblolly pine and are listed in Table S14 (8).

SI Note 3. Forest SOC modeling
RothC model uses five conceptual carbon pools to simulate SOC turnover, including four active pools (Resistant Plant Material (RPM), Decomposable Plant Material (DPM), Microbial Biomass (BIO), and Humified Organic Matter (HUM)) and one inactive pool (Inert Organic Matter (IOM)) ( Fig. S10) (2). As shown in Fig. S10, carbon from plant materials is divided into DPM and RPM based on a material-specific DPM/RMP ratio. Then CO2, BIO, and HUM are formed from the decomposition of DPM and RPM according to CO2/(BIO+HUM) ratio and BIO/HUM ratio. The CO2/(BIO+HUM) ratio relies on the soil clay content, while the BIO/HUM ratio is 0.85 (2). More CO2, BIO, and HUM are formed from BIO and HUM. The equations for calculating the amount of organic carbon in a carbon pool that decomposes in a particular month can be found in the model user document (2). RothC can be run with two modes: forward and inverse. For the forward modes, changes in SOC are calculated for equilibrium state or short-term with known monthly plant carbon input and other input data, including monthly precipitation, monthly open pan evaporation, monthly mean air temperature, clay content of soil, DMP/RPM ratio, soil cover state, monthly manure input (if any), and depth of soil layer. For the inverse modes, the model can generate monthly plant carbon input at an equilibrium state with known total SOC content and other input data discussed above (2).
This study simulated SOC changes for a century timescale with initial SOC contents of the five carbon pools and known monthly carbon input from tree materials that was calculated by the forest growth and yield model. However, the data of five carbon pools are rarely known from existing soil databases. To build the initial SOC contents for the five carbon pools, we first calculated the quantity of plant carbon monthly required by soil to maintain an equilibrium state using the inverse mode of RothC and present total SOC content data collected from existing soil databases. In year zero of the study timeframe, the SOC is assumed to be at an equilibrium state. The results of monthly plant carbon input generated from this step were then used to estimate the initial SOC contents for the five carbon pools in year zero using RothC in forward mode. In step three, these modelled initial SOC content of the five carbon pools were used to run a short-term run (100 years) in the forward mode to simulate changes in total SOC content. Varied monthly plant carbon input from forest litterfall, biomass from thinning, snags, harvesting residues, and post-harvest roots under different scenarios were calculated by the forest growth and yield models discussed in SI Note 2 and used as inputs to run the model in step three. A DPM/RPM ratio of 0.25 is used for woodland (2). Other assumptions include: 1) 30 cm topsoil depth, 2) all year-around soil covered by trees, 3) no manure input.
For protection AR (baseline), we implemented an average value of 0.23 Mg C/ha/yr for annual soil carbon stock increase for baselines (10). For the calculation of GHG emissions from soil and forest floor for baseline, detailed calculations can be found in SI note 1.

SI Note 4. Lumber Production
The first unit operation in saw mills is debarking (11). Then logs without bark are sawn into green lumber, slabs/chips, and wet sawdust. The yield of wet lumber, slabs/chips, and wet sawdust from logs was estimated based on the literature data (12,13) and documented in Table S11. In this study, all the slabs/chips were assumed to be sold to produce durable wood products (12). Bark and wet sawdust, as well as dry shavings/chips and sawdust from the planing (planing process is discussed as follows), were used for energy generation in the saw mill.
Wet lumber are dried at 90-120 C (dry bulb temperature) in a dry kiln operating to reach the targeted moisture content (14,15). Mill residues (bark, wet sawdust, dry planing shavings/chips and sawdust) are combusted to provide the energy for the dry kiln. When mill residues are not sufficient for the energy demand, natural gas is used; when energy is excessive, mill residues are used for power generation (12).
The energy demand for kiln drying was determined as the total heat demand in drying divided by the overall energy efficiency for energy generation and drying. The data were collected from the literature and shown in Table S11 (14,16,17). The total heat demand was determined based on the the desorption heat (MJ) needed for evaporating 1 kg water out of wood (18). The total energy demand for kiln drying is met by the total lower heating value (LHV) of mill residues, given by equation (S16) (19,20). HHV is higher heating value (assumed 20 MJ/oven dry kg (odkg) for wood and 20.5 MJ/odkg for bark) (12); MC is moisture content (dry basis) of mill residues; H is hydrogen content percentage in fuel (assuming 6) (16,(19)(20)(21). The LHV for natural gas is 47.1 MJ/kg according to (22).
Then the dried lumber is planed to remove uneven surface and generate finished lumber (23). Dry shavings/chips and sawdust generated in this process are viewed as mill residues (16,23). Then dimensional lumber is ready for distribution to the market.
As a structural material, the life span of lumber is assumed as 30 years based on the report by the National Association of Home Builders and Bank of America Home Equity (24). Then lumber is assumed to be landfilled. The calculation of end-of-life GHG emissions is shown in SI Note 5.

SI Note 5. Landfill of wood waste
Landfilled wood waste emits GHG emissions (majorly CO2 and CH4) through a decay process (25)(26)(27). Due to the high GWP factor of CH4 (27.0 for non-fossil CH4 GWP-100) and potential energy recovery value of CH4, landfill gas recovery for energy generation is in trend in the US according to the US Environmental Protection Agency (EPA) (28). In this study, the GHG emissions of landfill are estimated in two steps. First, CH4-rich GHG emissions from landfill decay were estimated based on the Intergovernmental Panel on Climate Change (IPCC) First Order Decay method (12,25). Second, part of the landfill gas is recovered and combusted for power generation.
All the parameters and values are collected from the literature and documented in Table S15. The GHG emissions from landfilled lumber and CLT are estimated using the same methods presented in this section.
Equations S17 and S18 show the method of the IPCC First Order Decay (25).
In equation S17, Cdecomposed is the accumulative decomposed carbon mass from year 0 to year t. W is the mass of deposited wood waste (wet basis); DOC is the degradable organic carbon of wood waste; DOCf is the faction of DOC that can decompose; k is the landfill decay rate (29). In equation S18, MCF is the CH4 correction factor which is determined by site management (e.g., disposal depth in the soil, anaerobic conditions) (25); F is the volume fraction of CH4 in landfill gas; R is the total recovered CH4 portion by energy recover device (29). This study adopts the typical recovery value of R to be 0.75 based on Anshassi et al. (30). OX is the average oxidation factor describing the fraction of oxidized methane (29). After the CH4 emission is determined, this study estimated the CO2 emissions of wood waste landfills using the experimental data on the volume rate of CH4 (before being recovered) to CO2 (31). The recovered landfill gas is combusted to generate power. The generated power is calculated based on the total LHV of landfill gas and electricity generation efficiency (30). Then the substitution credit for the generated power is considered to replace the US SERC region market mix using the data from ecoinvent 3.6 cut-off database (32). The detailed values are shown in Table S15.

SI Note 6. CLT production and potential substitution benefits
At the CLT plant, lumber preparation ensures the quality of CLT production, including grading, grouping, and moisture content screening (33)(34)(35). In lumber preparation, the boards undergo visual grading and grouping process. The preparation groups the lumber for different layers and directions, since the lumber in longitudinal layers need to reach visual grade No. 2 and the lumber in transverse layers need to reach visual grade No. 3 (36). Then the lumber moisture content is measured. In this study, the dimensional lumber from the saw mills has reached the target moisture of 12%, therefore additional drying is not necessary. Then the selected and grouped boards are longitudinally end-jointed to make long, continuous lumber for layering and gluing (40). Fingerjointing is assumed to be the end-jointing type (12,37). The longitudinally assembled lumber is then layered, glued with resin, and pressed to form CLT panels (38). The resin in this study was selected to be melamine formaldehyde (MF), a commonly adopted resin in finger-jointing and facebonding (39,40). After pressing, the finishing steps are planing and end cutting. Planing can remove excessive resin and finalize uneven surfaces. End cutting is finished by Computerized Numerical Control (CNC) to reach the customized shape (41). CLT is then packaged and ready for transport to the construction site. The waste produced in the CLT plant are landfilled. The data (e.g., energy and chemical usage) for this stage were collected from the literature and documented in Table S12.
This study includes the potential substitution benefits of CLT replacing traditional reinforced concrete and steel as building materials (42). The structural material consumption for CLT buildings and steel & concrete buildings were estimated based on the work by D'Amico et al. (42). Table S16 shows the average structural material usage for two types of buildings per m2 floor area (42,43). Then the processes from ecoinvent database (see Table S10) were used to determine the GHG emission reduction per m 2 floor area, and the results were normalized to 1 metric ton CLT panel used (1.13 t CO2e t -1 CLT panel used). For steel, after the usage, the end-of-life recycling rate is assumed to be 35% based on the ecoinvent 3.6 cut-off database (43). The remaining steel is landfilled. The upstream production burdens of steel and concrete products and landfill of steel were derived from ecoinvent 3.6 cut-off database (see Table S10) (32).

SI Note 7. Biochar Production
The flowchart of the biochar plant is shown in Fig. S7. In pretreatment, initial grinding reduces the feedstock size to ~50 mm (44). Then the feedstocks are dried in the rotary drum dryer to reach the moisture content 10% (green weight basis) in an ambient temperature range of 160 C to 180 C (45). The dried biomass is further ground in the hammer mill to reach 2.5-3.8 mm. Biomass is then pyrolyzed in a fluidized bed reactor at 500 C and 1 atmosphere for 60 minutes in the nitrogen ambient (46). Pyrolysis outputs include biochar and gaseous products. Biochar is separated from the gaseous products that are ducted to the combustor to provide heat for pretreatment and pyrolysis. Pyrolysis kinetics were simulated in Aspen Plus through the multistep reaction mechanism method (MSRM) (46,47). In the simulation, the biomass feedstock is first decomposed into major lignocellulosic components (cellulose, hemicellulose, and lignin) and then goes through a series of reactions (48). This study selects glucose for cellulose, xylose for hemicellulose, lignin-C, lignin-O, and lignin-H for lignin (48). As the biomass composition data were collected from the ultimate analysis (e.g., carbon content, oxygen content), the triangular method was applied to calculate the mass fraction of each model compound (49). In total, four reactors were used in the simulation. This first reactor (RYield) breaks the biomass into five model compounds and ash, the second reactor (RBatch) simulates the primary pyrolysis kinetic reactions, the third reactor (RCSTR) simulates the tar cracking reactions, the fourth reactor turns the remaining metaplastic components into biochar (50).
Then multi-stage cyclones separate the biochar from the gaseous product. The gaseous products from the pyrolyzer are combusted to produce heat for the rotary drum dryer and pyrolysis reactor. If the heat from the gaseous product is not sufficient, natural gas will be combusted. The pyrolysis reaction parameter settings are provided in Table S17. The electricity and diesel consumption of unit operations in the biochar plant and the other inputs are displayed in Table S17.

SI Note 8. End-of-life of biochar
This study used the three-pool exponential decay method presented by Woolf et al. (51) to model the GHG emissions from the slow decay of biochar after application, as shown in equation (S19). cremain is the remaining carbon mass of biochar after year t for the initial carbon input cinitial (51). Cinitial is derived from the Aspen Plus simulation. Three carbon pools in biochar are C1, C2, C3, respectively, with corresponding decay rate k1, k2, k3. The value of C1, C2, C3, k1, k2, k3 are derived from the regression data of experimental results (51,52) and displayed in Table S4. FT is the ratio of target temperature T to the reference temperature Tref and is a function of Q10 as shown in equation (S20). Tref in this study adopted the literature data of biochar decay experimental results (pine wood-derived biochar) (52) (see Table S4). T is the annual average temperature of the area where biochar is applied. As mentioned in the main text, this study assumes that derived biochar is applied in similar climate conditions with forest locations. Hence, T adopts the same annual average temperature of forest locations. Q10 describes the impacts of temperature on the decomposition rates of biochar carbon (51) and is derived by equation (S21).
For the whole area studied in this work, averagely, the mean residence time (MRT = 1/k) of the biochar is 550 years (53,54). However, the MRT of the topsoil SOC in the forests highly depends on climate conditions such as temperature and precipitation. MRT of SOC is generally lower in low-latitude zones and increasing toward high-latitude zones. One estimated range of mean MRT of soil organic carbon in the study area (N30°-N40°) is 14-16 years (55), another indicates a range of 12-29 years for temperate forests (56), and some others show a range of 10-30 years for the southeastern US (57,58), which are all much shorter than the MRT of biochar. If forest residues are left in soil to form SOC instead of being converted to biochar, it is expected that the biogenic carbon stored in SOC will resident much shorter time than the carbon stored in biochar.

SI Note 9. Uncertainty analysis
Uncertainty analysis of this study is addressed by the equations S22-S23 below.
and are the pessimistic and optimistic values of the results (i.e., net carbon stock and net GHG balances) that are impacted by two system, namely forest system and product system. and are the pessimistic and optimistic carbon benefits of the forest system; and are pessimistic and optimistic values of product-level gateto-grave (including production and end-of-life) GHG emissions.
and are derived from equations (S24) and (S25). and are the minimum and maximum values of forest carbon yield (Mg C/ha) that describes the carbon yield in harvested logs, residues, litterfall, and dead roots. is determined by equation (S26) as shown below.
describes the biomass yield from the forest and generated based on the methods described in SI Note 2.
is the carbon content (as a percentage) of loblolly pine that is estimated at the 95% confidence interval based on literature data (60-76) (see SI Tables S1 and S2).
In equations (S24) and (S25), and are the minimum and maximum value for the simulated SOC stocks from RothC, respectively. For and , the minimum and maximum SOC in uncertainty analysis is calculated by equations (S27)-(S28) derived from the method by the Food and Agriculture Organization (FAO) of the United Nations (59).
Where and are repectively the minimum and maximum value of the initial SOC stocks estimated at the 95% confidence interval based on the variation within the aggregated 1km×1km grid cells (considering the original values from 30m×30m grids in the ISRIC-World Soil Information); and are the sources of forest carbon inputs to the soil, including carbons from litterfall, biomass from thinning, snags, harvesting residues, and post-harvest roots; and are respectively the minimum and maximum value of monthly mean temperature estimated at the 95% confidence interval based on the interannual variation between year 1991-2020 from the CRU TS 4.05 dataset. and are respectively the minimum and maximum value of monthly precipitation estimated at the 95% confidence interval based on the interannual variation between year 1991-2020 from the CRU TS 4.05 dataset. and are respectively the minimum and maximum value of soil clay content estimated at the 95% confidence interval based on the variation within the aggregated 1km×1km grid cells (considering the original values from 30m×30m grids in the ISRIC-World Soil Information). For and , as stated in the main text, this study deploys the two-step uncertainty analysis for LCA for lumber, CLT, and biochar. The first step is to conduct the sensitivity analysis for the life-cycle GHG emissions for each individual product (i.e., lumber, CLT, and biochar) (the parameters are shown in SI Table S8, S11, S12, S15, S16, and S17). Then those parameters with impacts larger than 5% are selected. The ranges of these parameters used in LCA for lumber, CLT, and biochar, are collected from the literature and shown in SI Table S1. After the sensitivity analysis, there are four parameters selected, namely mill waste recovery rate ( ), landfill decay rate (k), landfill MCF (MCF), and potential material substitution for steel frame in traditional buildings by CLT ( ). The second step is to perform the uncertainty analysis for and based on equations (S29)-(S30).
After deriving , , , , the uncertainty analysis is performed based on equations (S22)-(S23).   (A) Low density without thinning for S1 -Lumber without forest residues removal, (B) Low density without thinning for S2 -CLT with 50% forest residues removal for biochar, (C) Low density without thinning for S3 -CLT with 100% forest residues removal for biochar, (D) Low density with thinning for S1 -Lumber without forest residues removal, (E) Low density with thinning for S2 -CLT with 50% forest residues removal for biochar, (F) Low density with thinning for S3 -CLT with 100% forest residues removal for biochar, (G) High density with thinning for S1 -Lumber without forest residues removal, (H) High density with thinning for S2 -CLT with 50% forest residues removal for biochar, and (I) High density with thinning for S3 -CLT with 100% forest residues removal for biochar. Note: the GHG emission of CLT EOL and Biochar decay are not visible due to their minimal values compared to others. These figures are based on average values.  (A) Low density without thinning for S1 -Lumber without forest residues removal, (B) Low density without thinning for S2 -CLT with 50% forest residues removal for biochar, (C) Low density without thinning for S3 -CLT with 100% forest residues removal for biochar, (D) Low density with thinning for S1 -Lumber without forest residues removal, (E) Low density with thinning for S2 -CLT with 50% forest residues removal for biochar, (F) Low density with thinning for S3 -CLT with 100% forest residues removal for biochar, (G) High density with thinning for S1 -Lumber without forest residues removal, (H) High density with thinning for S2 -CLT with 50% forest residues removal for biochar, and (I) High density with thinning for S3 -CLT with 100% forest residues removal for biochar. These figures are based on average values.      Tables   Table S1. Parameter ranges used in the sensitivity analysis.